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We used methods of Bayesian statistical inference and the principle of maximum entropy to an- 
alytically continue imaginary-time Green's function generated in quantum Monte Carlo simulations 
to obtain the real-time Green's functions. For test problems, we considered chains of harmonic and 
anharmonic oscillators whose properties we simulated by a hybrid path-integral quantum Monte 
Carlo method. From the imaginary-time displacement-displacement Green's function, we first ob- 
- - - tained its spectral density. For harmonic oscillators, we demonstrated the peaks of this function were 

I in the correct position and their area satisfied a sum rule. Additionally, as a function of wavenum- 

ber, the peak positions followed the correct dispersion relation. For a double-well oscillator, we 
demonstrated the peak location correctly predicted the tunnel splitting. Transforming the spectral 
densities to real-time Green's functions, we conclude that we can predict the real-time dynamics 
■ for length of times corresponding to 5 to 10 times the natural period of the model. The length of 

' time was limited by an overbroadening of the peaks in the spectral density caused by the simulation 

, algorithm. 
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; I. INTRODUCTION 

o 

, One of the goals for doing computer simulations is the production of information useful in the interpretation and 
T— I ■ design of experiments. Notwithstanding important issues regarding Hamiltonian selection and parameterization, the 
interface of simulations with experiment is particularly challenging for quantum systems. The current Monte Carlo 
algorithms, whether they impose quantum particle statistics constraints or not, are performed either in real-time t 
or in imaginary-time (Euclidean time) r — —it. In real-time, the propagator exp(zti/) for a system, described by a 
Hamiltonian H, oscillates wildly at long-times. Analytically, these rapid oscillations self-cancel, but a Monte Carlo 
process, as it is typically used, has difficulty achieving this cancellation. As a consequence, modifications of the basic 
algorithms have been proposed to extend the simulations as long as possible in the real-time domain With these 
^ ■ new algorithms, simulations typically produce dynamics extending to 2 to 3 times the natural periods of the systems. 
O , In imaginary-time, the propagator exp(— rTJ) is diffusive and the rapid oscillations are avoided. Correlations functions 
O • G(r), however, are now a function of imaginary-time, and such functions do not easily convey the actual dynamics 
of the system. In principle, real-time correlation (Green's) functions G{t) can be obtained from the imaginary-time 
ones by the process of analytic continuation. In practice, this process is difficult because it is ill-posed and because 
^ , the Monte Carlo data is incomplete and noisy. 

' Recently, procedures were proposed to perform this analytic continuation Q. These procedures draw heavily 
upon methods of Bayesian statistical inference and the principle of maximum entropy to infer from imaginary-time 
correlation functions their associated spectral densities A{u}). Through linear-response theory, the spectral densities 
represent the spectra associated with numerous real-time measurements of current-current, spin-spin, etc. correlations 
functions. What apparently has not yet been tried is performing the Hilbert transform of these spectral densities to 
obtain the frequency-dependent retarded correlation function and then Fourier transforming this quantity to obtain 
the real-time correlation function. In this paper, we will carry out these additional steps. By doing this, we hope to 
gain a greater understanding of the physical content present in the spectral density returned by the Bayesian methods. 
We expected that the resultant real-time information would be limited by the approximate and probabilistic nature 
of the analytic continuation methods. We found, however, that the distance in real-time over which our results are 
valid is limited primarily by the ability of the simulation algorithm to produce good data. To interface profitably with 
the numerical analytic continuation, the simulation algorithm has to produce high quality data consistent with the 
assumptions of procedures. The algorithm we used had problems doing this, and we will describe the measures taken 
to reduce this difficulty. Even so, in most cases we were able to extend in real-time up to factor of 10 natural periods 
of the physical systems. Longer extensions are possible and require longer Monte Carlo runs. For present purposes, 
we had no physical motivation to do so. 

In Section II, we will discuss the various models studied. We simulated a particle moving in single harmonic and 
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double- well anharmonic potential and a collection of particles moving in chains of these potentials. For these models 
we know the exact solutions. By calculating their properties numerically, we can benchmark our methods. Certain 
properties of a single double-well potential, like the tunnel splitting, are easily obtained numerically. The phase 
diagram for a chain of such oscillators is also known This type of chain can exist in a symmetric or displacive 
state. In Section III, we summarize the numerical analytic continuation procedure we used and discuss our simulation 
technique. Modifying the simulation technique to be more naturally ergodic and to produce data with short statistical 
auto-correlation times was the most difficult and restrictive part of our study. We present our results and conclusions 
in Sections IV and V. 



II. MODELS 



We simulated five Hamiltonians. One was that for a single harmonic oscillator 
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which has the natural frequency = vtT™- Another was that for a chain of N such oscillators 

i 

Fourier transforming the displacements Xi and momenta p^, we can of course rewrite this second Hamiltonian as 

fc 

where k = — tt, —{N — 1)-k/N, . . . , tt and uo'j. — 2loq{\ — cos fc). In this form, the Hamiltonian is explicitly expressed as 
a collection on N independent simple oscillators. The natural frequency of an oscillator is Wfc. The third Hamiltonian 
was a variant of the harmonic chain 

i 

which after Fourier transforming the displacements becomes 

fc 

where = 1 + 2ti)g(l — cosfc). In this form, the Hamiltonian is again explicitly expressed as a collection on N 
independent simple oscillators but with a dispersion relation that has a non-zero frequency at fc = 0. The two other 
Hamiltonians were a single, symmetric, double-well potential 

+ - If (6) 

which has well bottoms at a; = ±1 and a barrier height of unity at a; 0, and a chain of such potentials 

^ = # + ?(^' - ^-i)' + -M - 1)' (7) 

^ — ' 2m 2 4 

i 

In the chain, e sets the barrier height of every double-well. In both chains we assumed periodic boundary conditions 
and disallowed particle exchange. 

For these Hamiltonians, our simulations obtained estimates of the imaginary-time displacement-displacement 
Green's function 

G,,{t) ^ G,^,{t) ^ {TrX,{T)x,m (8) 
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Here, the angular brackets denote thermal averaging. It is more convenient and illuminating to work with the spatial 
Fourier transform Gkir) of and it is known Q that 

^ ' 27r 1 - e-/3" ^ ' 

where Ak{uj) is the spectral density. This function has the properties that 

Ak{io) = ~Ak{~uo) (10) 
The odd symmetry of Ak{uj) allows us to re-express @ as 

and it is straightforward to show that Ak{uj) obeys the sum rule 

The spectral density Ak{Lo) is also related to the frequency Fourier transform G^'{lo) of the real-time, retarded 
Green's function Q 

G^{t) ^ -zemx.k{t),xkm) (13) 



via 



27r 7-00 w — w' -|- iry 
where < 77 ^ 1, from which it follows that 

1 r°° 

G§{t) ^ ~i0it)— duAkiio)e-''^'e-^* (15) 

For an individual harmonic oscillator of frequency ujk, the eigenstates and energies are exactly known, and all the 
quantities in the above paragraph are known analytically: 



Gkir) = r-r-4 7K^C(^^^ 

Zmujk smh(pCi;fe/2) 



(16) 



Ak{uj) = Tr[6{uj - ojk) ~ 5{uj + ujk)]/mujk (17) 

G«(c.) = -J— ( 1 1 ^ (18) 

2mLJk \uj — uJk + iij U! + ojk + irj J 

Gj!it) = 0{t)—sm(LJkt)e-'^' (19) 
muj 

(4) = ^coth(/3c^,/2) (20) 

The eigenstates and energies of a single harmonic oscillator have definite, well-known characteristics: Because the 
potential is symmetric about x = 0, the eigenfunctions have alternating parity. The ground state has even parity and 
an energy ^cok- The energies of the excited states are regularly spaced at intervals of tUk- The double- well potential 
is also symmetric about x = 0, and its eigenstates also alternate in parity with the ground state again having even 
parity. The precise details about the energy spectrum, however, are only available numerically. When these states 
lie below the barrier, and particularly for deep wells, they group into widely-separated, nearly degenerate pairs. The 
separation in energies within and between pairs is called the tunnel spitting. The spectral density is dominated 
by terms with matrix elements involving states and 1. Contributions from matrix elements of x involving (0,3), 
(3,4), (0,5), (3,5), (5,6), etc. have smaller contribution to the spectral density because of smaller overlap between 
the eigenstates. Additionally, most can be "frozen-out" by making the temperature at least comparable to E3 — Eq. 
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This temperature range is the one in which we generally worked. The nature of the energy levels and eigenstates is 
schematically represented in Fig. 1. 

The spatial Fourier transformation of the Hamiltonian of the double- well chain does not produce as system of inde- 
pendent oscillators. This is the essence of its non-linearity. The model, however, has an interesting zero-temperature 
phase diagram as a function of the parameters e and 7 |^ . Roughly, e is a measure of the barrier height relative to the 
frequency of inter-site oscillation and the frequency of oscillation associated with the well-bottom. When the barrier 
height is large, the particles collectively are displaced to the left or the right of their classical equilibrium positions in 
a broken symmetry state characterized by a non-zero value of the mean-squared displacement. As the barrier height 
is lowered, a critical value is reached where because of zero-point motion and tunneling, the particles collectively 
make a transition to a state where the mean-squared displacement of each is zero. Accordingly, a simulation of the 
chain, done in one of these two thermodynamic phases, is expected to exhibit different quantum characteristics in the 
spectral density. 



III. METHODS 



A. Hybrid Path-Integral Monte Carlo Method 

Our Monte Carlo simulations will be based on the Feynman path integral formulation of quantum mechanics. In 
imaginary-time, this formulation represents the partition function Z as 



I?a;e-'5[^(*)1 (21) 



where 



S\x{t)\ = / dTH\x(T) 

Jo 



(22) 



is the action corresponding to the path x{t) and H[x{t)] represents the path dependence of the Hamiltonian. The 
Monte Carlo method is used to perform the integration over the paths in (|2^). It is applied after the integral in ( [2^ ) 
is approximated by a sum over L steps in imaginary-time each of length At and the momenta are approximated by 
a forward-difference approximation between successive displacements in imaginary-time: 

^x^{T) Xiir + Ar) - x^jr) 
P^{T) = "1^— « m — (23) 

For a one-dimensional system of N particles, the action becomes 

N.L 

S = ArY: ^ C'''^;'^' f + V{x.,,x.^,,) (24) 

ij = l 

This form is similar in appearance to a classical two-body potential energy defined on a A'' x L lattice where at a 
given r the particles interact through the potential energy function of the original problem (scaled by At), and at a 
given position they interact by a harmonic potential with a spring constant equal to m/ At. For a single particle, the 
summation over the spatial coordinate i is dropped and the discretized action has the interpretation of a chain where 
the particles at imaginary-time move in a potential ATV{xj) while interacting with particles at neighboring times by 
a harmonic potential with spring constant 1/Ar. For the models we are considering, the discretized actions are: 

1. Single Harmonic Oscillator: 
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Xj + i 



At 
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(25) 



2. Harmonic Chain: 



S = Ai 



N,L 

Em 



At 
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(26) 
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3. Harmonic Chain plus on-site oscillator: 



S = Ai 



N,L 

Em 



+ 2 



\2 . _|_ 2 



(27) 



4. Single Double- Well Potential: 
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(28) 



5. Double- Well Chain: 



S = eAi 



L.N 

Em 



At 



(a; 



(29) 



The simplest way to apply the Monte Carlo method is to move repeatedly from point to point on the space-time 
lattice, at each point propose a change in the coordinate, and accept the change via the Metropolis 



algorithm with probability min[l, exp(— AS)] where A^* is the change in the value of the action proposed by the 
proposed change Q . This method is often called path-integral Monte Carlo (PIMC) . 

An alternative to the Monte Carlo evaluation of the path- integral is a molecular dynamics evaluation Here, a 
fictitious momentum TTi., is associated with each point to define a pseudo Hamiltonian 



Hp — 



TT, 



N,L 

y 

^ 2m 



(30) 



and standard molecular dynamics techniques are used to sample phase space. The approach takes advantage of the 
classical nature of the fields in the path-integral formulation and produces the correct statistical mechanics because 
in classical statistical mechanics the momentum degrees of freedom can be integrated out of the partition function. 
The method is often called path-integral molecular dynamics. 

At the core of the method we used is the hybrid Monte Carlo approach suggested by Duane et al. 0] which 
combines the molecular dynamics approach with the Monte Carlo procedure to obtain the best features of both 
methods. The general expectation is faster equilibration of the simulation and shorter auto-correlation times between 
measured quantities. With this method, the following steps are cycled: For a given set of Xij the corresponding 
pseudo momenta are assigned values randomly from a Maxwell-Boltzmann distribution for the velocities. The energy 
is computed. Next, both the momenta and displacements are evolved by molecular dynamics for some pseudo time tp. 
The energy is recomputed. Then the evolved displacements are accepted with probability min[l, exp(— Ai^)], where 
AE is the difference in energy between the initial and final configuration. Normally, molecular dynamics is energy 
conserving so the evolved displacements would always be accepted. A guiding idea behind the hybrid method is to 
use for the molecular dynamics with emphasis on fast integration, as opposed to accurate integration, and to adjust 
the size of these steps At and tp so the Monte Carlo decision accepts 90 — 95% of the configuration. The molecular 
dynamics method globally updates all the displacements and is a computationally efficient procedure. The Monte 
Carlo procedure maintains detailed balance to ensure proper equilibrium averages and filters out the results of "bad" 
integrations. 

We found, however, that this simple form of the hybrid method was inadequate for present purposes. The output 
of our simulation is to be used as input to maximum entropy procedures to execute the analytic continuation. As 
we will discuss below, the analytic continuation problem is an ill-posed problem and hence is very sensitive to the 
size of the errors associated with the input data. For a fixed amount of computer time, reducing the size of the error 
efficiently by a Monte Carlo method requires shortening the auto-correlation times between measurements. In our 
computed Green's function, in spite of small estimates for the error bars, we would often see small (within the error 
bars) unphysical sawtooth-like structure in regions about r — (3/2. Following a simple procedure suggested by Neal 
P, we could generally remove this structure and also be more ergodic. His suggestion was after each Monte Carlo 
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decision to reverse the direction of the molecular time integration, i.e., At —At, with probability 1/2. A smaller 
improvement is achieved by not using a fixed length for the time integration in the molecular dynamics simulation but 
rather to choose the length randomly from the interval {tp — 5,tp + S) where tp and S are chosen so the Monte Carlo 
acceptance rate is in the 90 to 95 per cent range. We will refer to the combined method as the hybrid path-integral 
Monte Carlo method (HPIMC). Adjusting the Monte Carlo acceptance rate is not the entire story. First, it seems 
best to insure tp is several times larger than the natural period associated with the slowest significant modes in the 
systems and then choose S to fix the acceptance ratio. 

We remark that Fahy and Hamann ||^ observed for the standard hybrid method the existence of a critical time tc 
(dependent of model parameters) demarking non-ergodic and chaotic behavior in the results of the time integration. 
For a harmonic system, tc is infinite which suggests the inapplicability of the method to a harmonic system. We 
observed the behavior they found but whether tp was larger or smaller than t^ had only small consequences on our 
measured results. As we illustrate below, we achieved very accurate results for the harmonic models. 

We also remark that the results of the simulations depend on the size of At. By performing simulations for several 
different values of At, we could in principle extrapolate the results to the At = limit. We did not do this but 
instead observed that simulations performed with different values of At gave very similar results. 



The maximum entropy method | 
we rewrite this equations as 



where the kernel 



B. Maximum Entropy Method 

is used to regularize the solution of (pi). Dropping the subscript k for convenience. 



G(t) = / dLuK{T,uj)[A{uj)/uj] 



Kir) 



1 w[e-^^ + ( 
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(31) 



(32) 



Because both the kernel and A{u!)/lu are regular a,t u — 0, we solve for A{u!)/llj and then trivially find A{llj). For 
discrete values of t and ui we approximate (^l]) as 



(33) 



^,3 



where Gi — G{Ti), Kij — K(Ti,ujj), and Aj = A{ujj)Aujj /u 



The Monte Carlo method will return estimates Gi of Gi and estimates of the sample variance cr^ on the Gi . With 



this information, a natural solution path to Aj would be to find the values of Aj that minimize 



E 



G, - ^.j^, 



(34) 
(35) 



This approach, however, almost always fails. One reason is that it ignores the strong correlations that normally exist 
between the measured values of Gi, i.e., values of the Green's function at different imaginary-times. At the very least, 
we must modify (B5|) to be 



X 



^(Gi - Gi)[C ^]i,jiGj - Gj 



(36) 



where Cij is the measured covariance among the values of Gi. The i-th diagonal element of G is simply af. This 
modification of the definition of while necessary is insufficient. The difficulty is the inverse problem, that is, solving 
( ^ ) for A{uj), is ill-posed. This condition is caused by the exponential character of the kernel at large values of co. 
At large large variations in A{u;) make little change in G(t). The simulation, on other other hand, gives noisy 
and incomplete information about G{t), and hence for a given set of Gi, an infinite number of Aj will satisfy the 
least-squared estimate (|36|). 
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The next level of solution seeks to regularize the minimization of by constraining it, i.e., minimizing 



x2-5]a,/,(A) (37) 

i 

where the are Lagrange multipliers and the fi{A) are functions of Aj representing possible constraints on the 
solution. Typical constraints include smoothness, non-negativity, sum rules, moments, etc. The difficulty with this 
approach is choosing the Lagrange multipliers. Ad hoc choices are commonplace. Often small changes in the values 
of these parameters produce massive changes in the results. 

The maximum entropy approach follows from the observation that the spectral density is interpretable as a proba- 
bility density function. The principle of maximum entropy states the probabilities should be assigned in such a way 
as to maximize 

S = Y,A,- m, - A, \n{A, /m, ) (38) 



Here, the rrij, called the default model, set the location of the maximum of S and the value of S at this point to 
le c 

entropy maximizes 



be zero. The default model is solution for Aj in the absence of other constraints on Aj. The method of maximum 



Q{A) ^aS^ (39) 

To fix a, an ad hoc procedure called historic maximum entropy is often used |lo|Jll[] . An more modern alternative 
is the Bayesian-based classic maximum entropy which uniquely determines a provided certain conditions are meet 
1^,0. Under these conditions solution for Aj is the most probable one. Unfortunately, these conditions seem often 
violated in the analytic condition problem. Accordingly, to estimate the Aj, we adopted a procedure suggested by 
Bryan. 

In Bryan's method p^ , for a given value of a, we find the A{a) that maximizes Q{A). For the solution to (p3|), we 
take 

A^ I da A{a)¥v[a\G] (40) 



where Pr[a|G] is the probability of a given the data G. Bayesian analysis shows that 

Pr[a|G] = PrH j VA ^^^^ (41) 

where Zj^ is the normalization factor for e~\^ , Zgip) is the normalization factor for e"'^, and Pr[Q!] is Jeffreys' prior. 
Details on the computation of this joint probability are discussed elsewhere [^. With this function, the integral (|40| ) 
is performed numerically. 

The most difficult part of the problem is not evaluating the maximum entropy equations but satisfying the statistical 
assumptions on which they are based. The principal assumption is 

Pr[G|A] = ^ (42) 

The meaning of this assumption is the measured values of Gi are statistically independent and distributed according 
to a multi-variable Gaussian distribution function defined by the covariance matrix C . Proper estimation of C is 
paramount. Under normal circumstance the data produced by the simulations do not satisfy these assumptions. The 
procedures we use to have the data approximate these assumptions are discussed elsewhere |^ . When we have proper 
data, our solution (^0|) usually shows good insensitivity to the choice of the default model. Additionally, the historic 
and classic maximum entropy solutions usually agree well with it. 

IV. RESULTS 

To determine spectral properties of models listed in Section II, we performed HPIMC simulations with up to 800 
bins of data (Nun) , each with up to 4000 measurements (Ngweep)- Simulations with large bin-sizes were necessary 
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to avoid nonergodic behavior of the HPIMC method when used for chains with double well potentials close to the 
zero-temperature phase transition point. Furthermore, we set the value for the imaginary time step to At = 0.25. 
This choice on the one hand was small enough to avoid errors associated by the discretization of otherwise continuous 
imaginary time scale r, and on the other hand was large enough to avoid unphysical correlations between successive 
imaginary-time measurements of the Green's function G(t). Since our calculations were performed at the inverse 
temperatures (3 between 1 and 10, the corresponding number of imaginary time steps L — /3At was between 40 and 
60. 

For a successful application of the HPIMC it is crucial to choose the proper value of the pseudo-time tp and the 
size of its step At in the molecular dynamic part of the simulation. Following Fahy and Hamman, we determined the 
value of the critical value for each case under consideration and then took tp > tc to avoid running the simulation 
in a non-ergodic regime. Typical values for tp were between 5 and 15 for double-well cases listed below. We stress 
that there was not much difference in the quality of the HPIMC data if we chose tp > tc or chose tp = tc/2. In 
addition, we obtained good data for the harmonic wells by choosing tp ~ 1/wo even though for this particular case 
tc — oo- If we define as the smallest nonzero frequency of the system, optimal values of the step size At are 
between 0.05/wo and 0.1/wo- Larger values of At lead to larger errors in the pseudo-time propagation which then 
lead to small acceptance ratios. Smaller values of At lead to longer compution times. Unless specified otherwise, we 
always chose m = =7 = 1. We emphasize that the HPIMC method is insensitive of the choice of the mass TOjt 
associated with the fictitious momentum nij. 

In Fig. 2, we show the displacement-displacement Green's function G{t), for a single harmonic oscillator, obtained 
by evaluating ( p^ for various values oiuoP, as a function of the imaginary-time variable r. These curves look similar 
to the Green's functions that we obtained numerically for this and the other models: For some parameters and 
temperatures, the G(t) varies little as a function of t; for others, it varies rapidly at the ends of the interval [0, (3) 
and is flat in the middle with values nearly equal to zero; and for still other parameters, it has a featureless, parabolic 
looking shape. 

The common features of these curves have several significant implications for the analytic continuation problem. 
First, we remark that from the quantum Monte Carlo simulations we obtain estimates of G(r) only at a relatively 
small number of discrete values r.;. The smoothness of the curves implies that the G{Ti) at neighboring values of 
are correlated. The computation of the covariance matrix in (|3^ ) is thus an important part of the of the analysis 
of the data. While the correlations among the different r values of G{t) make the interpretation of the assignment 
of an "error-bar" to a given r value delicate, such an assignment illustrates several difficulties inherent in the data 
that help to make the analytic continuation of the data often very difficult. In the case where G(r) is nearly flat, the 
errors bars mean that a number of values of G(r) are "within the error-bars" of each other. This situation, along 
with the correlations implied by sizable off-diagonal elements of the covariance matrix, means that only a subset, and 
often a small subset, of the measured values of G(r) represent independent data useful for the analytic continuation 
procedure. The analytic continuation near the classical limit can be very difficult. 

The situation with the rapid end-point variation and the flat nearly zero values is another difficult case. Again the 
fiat region generates a lose of useful values of G(r) and the smallness of G(t) in this region can engender situations 
where the error bars would imply that during the simulation estimates of G(t), which must be non-negative, were 
derived from ensemble values that included negative ones. The Monte Carlo algorithms in fact do not produce negative 
values but do produce highly skewed fluctuations about the mean. The Gaussian assumption for the likelihood function 
in ( ^2|) thus can only be approached in the limits of a large number of independent measurements when the central 
limit controls the data distribution. The ratio of the mean value to the estimated variance (signal to noise ratio) also 
indicates that the most effective data comes from those in the rapidly decreasing region. At low temperatures, the 
analytic continuation problem can become very difficult. 

The details of the simulation algorithm can also impact the quality of the results and data. In Fig. 3, we show 
G(t)/G(0) as a function of r obtained by two closely related simulation techniques for a single harmonic oscillator. 
We remark that the scale of abscissa is 1/100 of that of Fig. 2 and the ordinate shows r only in a narrow region at the 
symmetry point (3/2. The dashed curve is the analytic result obtained from ([l^). The data is represented by square 
markers were obtained by the HPIMC method without the time-reversed step; the data represented by the circles 
were obtained with the HPIMC method with the time-reversed step. In each case, the same number of Monte Carlo 
steps were made. One sees that the fluctuations with the HPIMC method without the time-reversed step are larger 
and that the error bars associated with the results suggest a dip into non-negative values of G(r). More signiflcantly, 
the results deviate from the exact curve by more than one standard deviation in the immediate vicinity of r = (3/2. 

The analytic continuation result for A{u}) from the data partially shown in Fig. 3 is shown in Fig. 4. From (pT|), the 
spectral density should be Q.b5{uj — 1). The solid curve is obtained from the HPIMC algorithm and shows a broadened 
delta-function at the right location with nearly the correct weight. The fraction of a per-cent difference from the 
correct weight is most likely a consequence of the small error caused by discretizing the imaginary-time derivatives. 
On the other hand, the dashed curve, which is obtained from the data obtained from the HIPMC algorithm without 
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the time-reversed step, is broader, located incorrectly, and has a larger discrepancy in its weight. The increased 
breadth is a consequence of the larger variance in the measured data. The incorrect location and poorer weight is a 
consequence of the small deviation from the exact value near r = /3/2, which in turn is a consequence of the Monte 
Carlo algorithm performing badly. 

In Fig. 5, with the data obtained from the HPIMC algorithm, we show the full analytic continuation form G(t) to 
G^{t). In Fig. 5a, the data (open circle) is compared with the exact results for G(t) (solid line) obtained from (|lq). 
In Fig. 5b, the dashed curve is the A{uj) obtained by the analytic continuation procedure, while the solid line is a 
Lorentzian at the same location. The real part of G{uj) is shown in Fig. 5c, where the dashed line is the quantum 
Monte Carlo results and the solid line is an analytic result obtained from the Lorentzian from Fig. 5c. The width 77 
of the Lorentzian shown in Fig. 5b was adjusted so the w = values of the two curves agreed. The single adjustment 
produced remarkably good agreement at high frequencies. The principal differences between the two curves are at 
uj — ±1 where divergences should exist as indicated by (|T^). Finally, G^{t) is shown in Fig. 5d. The solid line was 
obtained analytically and used the Lorentzian of Fig. 5b, while the dashed line is the quantum Monte Carlo result. 
The agreement between the results is satisfying and comparable in quality to that possible by real-time quantum 
Monte Carlo methods. The width of the Lorentzian was 7/ujo, i-c, about 7 times the natural period of the harmonic 
oscillator. This width is still controlled by the size of the variance in the measured values of G(t). 

For the case of the harmonic chain, we computed the Green's function Gfe(r) for each independent wave number 
and from it found the corresponding spectral density Ak{uj). As with the single oscillator, this function should be 
a delta-function located at and have a weight equal to Gfe(O) = (xl). As in the single oscillator case, instead of 
finding a delta- function, we found a Lorcntzian-like peak at the correct location with the correct weight. The peaks 
for different values of fc, however, had different widths. In general, the peak widths increased with increasing fc, 
correlating with the increased variance associated with the G^ (r) . The spectral densities for the three lowest values of 
k are shown in Fig. 6. The peak positions as a function of w give the phonon dispersion relation. Our determination 
of this relation is compared to the exact result in Fig. 7. The agreement is excellent. There is a difficulty that must 
be mentioned. At fc = and uik = 0, Gkir) is flat and the weight of the peak approaches infinity. For fc = 0, not 
surprisingly, we were unable to do the analytic continuation. This situation is why we considered the model defined 
by (^. Here, at fc = is not zero and the determination of Ak{uj), and subsequently tok, is possible for all values 
of k. The dispersion relation found from the analytic continuation agrees very well with exact result. 

For a single, double-well potential, the spectral density can give direct information about tunneling processes. 
This situation is illustrated in Fig. 8 where the spectral densities for the Hamiltonian, described by (|^), are shown 
for several different values of the parameter 7. It is straightforward to discretize Schroedinger's time-independent 
differential equation for this potential and find the eigenvalues of the resulting eigenvalue equation. We adjusted the 
model parameters so only a very few (usually one) of the lowest eigenstate lied below the barrier height, similar to the 
situation depicted in Fig. la. (For deeper wells, our simulation methods would get stuck in one well or the other for 
large numbers of Monte Carlo steps, and hence the algorithm effectively lost ergodicity.) In the cases reported, the 
temperatures of the simulations were also less than the separation between the two lowest lying pairs of eigenstates. 
Thus, our spectral densities only exhibited the transition between these two states, and the position of the peak 
gives a direct measure of the lowest frequency tunnel splitting. This position agreed very well with the exact value 
calculated from Schroedinger's equation. Additionally, the weight of the peak also agreed well with the sum rule (|l^). 

For a chain of double- well potentials (^, still a different situation presents itself. At zero temperature, depending on 
the values of e and 7, the system exists in either a broken symmetry phase, in which the mean-squared displacement is 
non-zero, or in a symmetric phase, in which the mean-squared displacement is zero. The model has a quantum phase 
transition. We found that simulation, if performed in the broken symmetry phase, stuck in one well or the other for 
large numbers of Monte Carlo steps, making it difficult to collect the amount of statistically independent information 
to do the analytic continuation. In Fig. 9, we report the Ak{uj) for two simulations done in the symmetric phase. The 
one in Fig. 9b is close to the zero-temperature phase boundary. As one moves closer to the T = phase boundary by 
increasing e, the k = peak moves towards uj = 0. This movement is a consequence of the decreasing probability for 
tunneling between the two minima of the double-well potential as the barrier height is increased. 

V. CONCLUDING REMARKS 

We used methods of Bayesian statistical inference and the principle of maximum entropy to analytically continue 
imaginary-time Green's function generated in quantum Monte Carlo simulations to obtain the real-time displacement- 
displacement Green's functions. For test problems, we considered chains of harmonic and anharmonic oscillators whose 
properties we simulated by a hybrid path- integral quantum Monte Carlo method (HPIMC). From the imaginary- 
time, displacement-displacement Green's function, we first obtained its spectral density. For harmonic oscillators. 
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we demonstrated that the peaks of this function were in the correct position and their area satisfied a sum rule. 
Additionally, as a function of wavenumber, the peak positions followed the correct dispersion relation. For a double- 
well oscillator, we demonstrated the peak location correctly predicted the tunnel splitting. Transforming the spectral 
densities to real-time Green's functions, we conclude that we can predict the real-time dynamics for length of times 
corresponding to 5 to 10 times the natural period of the model. The length of time was limited by the simulation 
algorithm and not directly by the analytic continuation procedure. 

The simulation algorithm influences the results in at least two ways. One way is the ease with which statistically 
independent, Gaussian-distributed data is obtained. In our experience, of the various quantum Monte Carlo methods 
we have used, path-integral Monte Carlo methods (PIMC) tend to produce data with long-ranged correlations, thus 
making statistical independence sometimes difficult to achieve. Achieving statistical independence is important for 
proper error estimation because of the sensitivity of ill-posed problems to errors in the data. With out statistical 
indepence these error are usually underestimated. That the data is Gaussian-distributed is necessary to satisfy the 
assumptions of the analytic continuation procedure. Currently, large amounts of data are used to force the simulation 
data to have proper statistical properties. The method of binned averages with many measurements in a bin, is 
used to achieve statistical independence. The central limit theorem is used to obtain a Gaussian distribution of the 
binned averages. Because of the low computational intensity necessary for the test cases considers, producing the 
necessary large amounts of data was not problematic. 

A second way the algorithm can influence the results is through broken ergodicity and unphysical results. Here, we 
are referring to the ragged structure we saw in G(r) with the hybrid method and to the small systematic difference 
between the exact and computed results. We have not previously seen similar problems. Small modifications in 
the simulation algorithm removed the problems but it was at first difficult to determine the source of the difficulty. 
Because the purpose of the research was not algorithmic development, we did not do any comparisons of HPIMC 
and other possible methods to determine relative efficiency and other merits. Again, the computation times for these 
simulations is small (a few hours on a modest workstation) and so we were unmotivated to make such comparisons. 
Other recently suggested approaches, e.g., |l^,|l^, should be considered as part of further studies. 

In general, we believe we have demonstrated that analytically continuing imaginary-time correlation functions, 
obtained from a quantum Monte Carlo simulation, to obtain real-time correlation functions is a feasible alternative to 
obtaining such real-time functions directly from a quantum Monte Carlo simulation done in real-time. One advantage 
in choosing this approach appears to be the longer length of time over which the real-time information is faithful to 
the correct result. Of course, this conclusion is based on one study of simple models; however, this study does strongly 
indicate the direction for further and more extensive work. 
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FIG. 1. Schematic representation of the energy levels in a harmonic and double- well potential. 



FIG. 2. Green's function G{t) for a particle in a harmonic-well, obtained from the analytical form (21) potential, as a 
function of imaginary-time t and at different values of /Jwo- 

FIG. 3. Comparison of analytical results for the Green's functions for a particle in an harmonic well (dashed curve) with 
numerical results (open circles and squares) obtained using HPIMC method (a) with and (b) without the time-reversal step in 
the molecular dynamic part of the algorithm. 

FIG. 4. Comparison of the spectral functions obtained from the data shown in Fig. 3. As in Fig. 3, in we present spectral 
functions extracted from HPIMC data (a) with and (b) without the time- reversal step. We also present G(0), the computed 
area under the curves. 

FIG. 5. Comparison of the Green's function (a) and its spectral fruntions for a particle in a harmonic well obtained by 

randomly changing the direction of the molecular dynamics time integration (open circles in (a) and full line in (b)) and by 
not randomly changing he integration time direction (open circles in (a) and full line in (b)). Presented is only a small portion 
of G{t) around t = P/2 where deviations from the analytic solution (dashed line in (a)) are enhanced. Also presented are the 
sum rules. Analytically, G(r = 0) = (x^) = 0.500. 

FIG. 6. Spectral functions for the harmonic chain at different values of k. Arrows above spectral functions denote exact 
positions of peaks. We also compare the numerically calculated sum-rules / with analytical values. 

FIG. 7. Dispersion relation ujk for harmonic chain with m = 1 and 7 = 1. Open circles correspond to discrete values of uJk 
for a 10-site system. 

FIG. 8. Spectral functions for a single double- well potential at different values of 7 and different temperatures 1//3. Arrows 
above spectral functions denote exact positions cuo of peaks as obtained by numerical solution of the double-well potential. 
We also compare sum-rules, calculated by integrating the spectral functions A{uj) over oj, with G(t = 0) =< >, obtained 
directly from the QMC calculation, u is the ratio 1/ijto where 77 was obtained by fitting A{u!) to the analytical form G(u!), see 
Eq. (23). 

FIG. 9. Spectral functions A{u}) for a chain of 20 sites at two different values of e: (a) e = 0.5 and (b) e = 0.75. We present 
the spectral functions for different wavevectors k in units of 27r/20. 
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